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ABSTRACT 

A new analytical model for constraining the extent of gravitationally bound struc¬ 
ture in the Universe is presented. This model is based on a simple modification of 
the spherical collapse model (SCM), and its performance in predicting the limits of 
bound structure in (V-body simulations is compared to that of two previous models 
with the aid of new software named COLDGaS - compute unified device architecture 
(CUDA) object location determination in GADGET2 snapshots - which was developed 
by the author. All of these models can be distilled down to a single unique parameter 
£, here named the critical parameter, which was found to have values of 3 and 1.18 
from the previous studies, and a value of 1.89 from the modified SCM. While still on 
the conservative side, this new model tends to better identify what structure is gravi¬ 
tationally bound in simulations. All three analytical models are applied to the Corona 
Borealis supercluster, with the modified SCM and £ = 1.18 model making predictions 
that are in agreement with recent work showing that A2056, A2061, A2065, A2067, 
and A2089 comprise a gravitationally bound supercluster. As an additional test, the 
modified SCM is used to estimate the mass within the turn around radius of the Virgo 
cluster, providing results in good agreement with studies relating the virial mass of 
clusters to the total mass within turn around. 
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1 INTRODUCTION 

With the release of the Planck survey data 
(Planck Collaboration et al. 2014), we have entered 
an unprecedented era of precision cosmology. The values 
of the cosmological parameters from that survey are in 
excellent agreement with the standard ACDM cosmology, 
meaning that it is all but certain that we live in a universe 
that has entered a phase of accelerated expansion. This 
will cause gravitationally bound structures, such as groups 
and clusters of galaxies as well as dense superclusters, to 
one day become ‘island universes’ (Nagamine & Loeb 2003; 
Araya-Melo et al. 2009), with other structures eventually 
expanding beyond the observable horizon. 

Locating the largest gravitationally bound structures 
has implications beyond simply defining objects des¬ 
tined to become ‘island universes’. Accurately determin¬ 
ing the number of extreme density peaks could place con¬ 
straints on the nature of the primordial fluctuation field 
(Sheth & Diaferio 2011; Park et al. 2012). So far, there 
are at least three extreme density peaks known within 
z ~ 0.2, the Shapley supercluster (SSC; Bardelli et al. 1993; 
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Proust et al. 2006), the Sloan Great Wall (Luparello et al. 
2011), and the Corona Borealis supercluster (CSC; 
Small et al. 1998; Luparello et al. 2011; Batiste & Batuski 
2013) all with estimated masses of ~ 10 16 /i -1 Mg. Recently 
Pearson, Batiste & Batuski (2014, hereafter PB") provided 
the most conclusive evidence to date that the CSC con¬ 
tains an extended gravitationally bound core consisting 
of five rich Abell clusters, A2056, A2061, A2065, A2067 
and A2089. Given the potential importance of these ob¬ 
jects, having a number of tools to help identify them is 
paramount. 

The spherical collapse model (SCM) has often been 
used to study non-linear evolution in the presence of a cos¬ 
mological constant (Lahav et al. 1991; Wang & Steinhardt 
1998; Nagamine & Loeb 2003; Percival 2005; Hoffman et al. 
2007). Of particular interest are two previous studies that 
attempted to define the limits of gravitationally bound 
structure in a universe undergoing accelerated expansion. 
Busha et al. (2003, hereafter B03) used the SCM, along 
with the assumption that galaxies would currently be re¬ 
ceding from overdensities with an unperturbed Hubble flow 
velocity, to develop an analytical model for the limits of 
bound structures in a ACDM universe. However, structures 
that are destined to be gravitationally bound will have likely 
slowed from the overall Hubble expansion at present, mean- 
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ing this model is quite likely to underestimate the extent of 
bound structure. 

Diinner et al. (2006, hereafter D06), seeking to improve 
upon the model of B03, applied a pure SCM to determine 
the limits of bound structure. Beginning with the Tolman- 
Bondi equation, and transforming it into dimensionless vari¬ 
ables, D06 were able to determine the radius of the ‘crit¬ 
ically bound shell’ without assuming anything about the 
present velocities of galaxies. This model was suggested as 
being more realistic since it should account for the slowing 
of ‘shells’ that are part of a bound structure. However, the 
pure SCM ignores the effects of external attractors which 
will be present in the real Universe. Since these external 
attractors will pull on the extremities of a structure, those 
parts may not actually be bound. This leads the model of 
D06 to be optimistic in its limits to gravitationally bound 
structure. 

By defining a parameter £, here named the critical pa¬ 
rameter, the results of both models can be expressed as, 


M ob j 
10 12 M 0 


^ £ ft 70 


\ 3 

■ 

1 Mpc J 


(1) 


Here, A/ 0 bj is the mass contained within the radius 
ro, and /170 = Ho/70 kms -1 Mpc -1 , instead of the usual 
h = Hq/100 kms -1 Mpc -1 . Thus, locating the limits of 
bound structure in the spherical approximation becomes 
the pursuit of that single parameter. B03 found a value 
of £bo 3 = 3, which has a certain appeal due to its simplic¬ 
ity, while D06 found a value of £do6 = 1.18. In an OCDM 
universe (A = 0), so long as the total energy of a shell, ki¬ 
netic plus potential, is negative, the shell will remain bound 
(Nagamine & Loeb 2003). However, to give some idea as to 
the sensitivity of £ to the value of A we can consider a shell 
which has a Hubble expansion velocity at present, i.e. the 
B03 model. Hoffman et al. (2007) found that for such a shell 
to have a turn around in an OCDM universe, the enclosed 
overdensity must be 5 C = 2.33, which corresponds to a value 
of £ = 0.57. 

In this paper, a more robust estimate of £ is sought us¬ 
ing different methods. First, the results of the IV-body sim¬ 
ulations performed by Pearson & Batuski (2013, hereafter 
P13) are used to obtain an estimate for £. Then, a simple 
modification to the standard SCM is proposed based on an 
observation of radial and tangential velocities in cosmolog¬ 
ical simulations. From that modified SCM, an estimate of 
£ is made in the same manner as in D06. The performance 
of these critical parameters is then tested with cosmologi¬ 
cal simulations performed with GADGET2, and with appli¬ 
cations to the CSC and the Virgo cluster. 

This paper is structured as follows: Section 2 presents 
the P13 simulation informed model. Section 3 presents the 
modified SCM. Section 4 is the comparison with the GAD- 
GET2 simulations. Section 5 demonstrates ways the analyti¬ 
cal models may be applied to real structures in our Universe, 
so long as redsliift independent distances are known. Lastly, 
in Section 6 the implications of the results are discussed. 

Throughout this paper it is assumed that f2 m ,o = 0.3, 
Da,o = 0.7, and h = Hq/100 kms -1 Mpc -1 . 


Table 1. CSC Critical Parameter Constraints. Column 1 lists 
the cluster pair being examined. Column 2 lists the mass of the 
larger cluster in the pair. Column 3 lists the initial separation, 
and Column 4 is the value of the critical parameter. 


Cluster Pair M n() j ro £ 

(10 15 h -1 M 0 ) (h- * 1 kpc) 


A2061/A2067 1.708“ 


A2065/A2089 0.951“ 


6326.80 

2.96 

8775.97 

1.11 

8417.99 

1.26 

6548.20 

2.67 

7426.02 

1.83 

9941.05 

0.76 

7039.69 

2.15 

6162.42 

3.20 

7928.76 

1.50 

7449.05 

1.81 

8173.69 

1.37 

8102.29 

1.41 

5690.67 

2.26 

5044.48 

3.25 


“Mass is listed once for reference. All values for a given pair use 
the same mass to determine £. 

2 P13 SIMULATION INFORMED LIMITS 

The simulations of P13 imply that the model of B03 would 
act as a lower limit to the true extent of bound structure, 
while the model of DOG would act as an upper limit. Given 
the excellent agreement of the SSC simulation results with 
the results of Munoz & Loeb (2008) as to the extent of 
bound structure, it is reasonable to assume that the sim¬ 
ulation results are reliable. 

The simulations of the CSC performed by P13 showed 
very little chance of there being extended gravitationally 
bound structure, due to using spectroscopic redshifts as 
distance indicators when the clusters had significant pecu¬ 
liar motions along the line-of-sight (Batiste & Batuski 2013; 
PB 2 ). However, there were some cases in which two differ¬ 
ent pairs of clusters showed up as bound, A2061/A2067 and 
A2065/A2089. It seems safe to then assume that in each of 
these cases, they must have been near the limits of bound 
structure. For this reason, they were examined to determine 
a value of £. Before that was done, the further requirement 
that the current kinetic energy of the pair be at least 95 per 
cent of the potential energy was applied, ensuring that the 
pairs were loosely bound. 

By assuming that the cluster pairs were at the limits 
of bound structure, the equality in equation ( 1 ) was taken 
and then solved for the critical parameter, 

t M °bj f 1 Mpc \ 3 

5 h 2 o (lO 12 M 0 ) \ ro J ■ U 

The initial separations and mass of the larger of the two 
clusters can be entered into equation ( 2 ), giving a value for 
the critical parameter. The results of applying this proce¬ 
dure to the CSC cluster pairs are shown in Table 1. Look¬ 
ing at the values of the critical parameters, we can see 
that there are two cases in which the criterion is more 
conservative than that of B03, meaning those are cases 
where the structure should definitely be bound. There are 
also two cases where the criterion is lower than that of 
D06, implying the clusters must have had motions towards 
each other, leading them to be bound. The rest lie some- 
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where between those two models. Taking a bi-weight me¬ 
dian (Beers, Flynn & Gebhardt 1990), a value of £ = 1.90 
is obtained. Confidence limits are determined using jack¬ 
knife resampling (Beers, Flynn & Gebhardt 1990), giving in 
the end £ S i m = 1.90 ± 0.37 with 95 per cent confidence. This 
confirms more rigorously the findings of P13, that the true 
limits to bound structure, at least in those simulations, lie 
between the B03 and D06 models. 


3 A MODIFIED SCM 


When it comes to applying the SCM to structures in the 
Universe, aside from the implicit assumption of spherical 
symmetry, there are two key weaknesses. First, the SCM 
excludes the effects of what can be called external attrac¬ 
tors, mass concentrations not bound to the central dominant 
cluster, which can pull galaxies away from the structure. 
Second, the SCM assumes that all motions are purely ra¬ 
dial. A quick examination of particles from simulation shows 
that there will be some fairly substantial tangential motions 
at present around structures that are still in the formation 
process. Fig. 1 shows the radial and tangential motions of 
particles around Object 1 from Simulation 2 (see section 4) 
at the present time (a = 1). The bottom panel shows that 
there are significant tangential motions all the way out to 
the turn around radius (~ ll/t -1 Mpc). Comparing with 
the top panel, we can see that it is reasonable to assume 
equal radial and tangential motions. 

In the standard SCM, assuming only radial motions, 
the energy equation is (Peebles 1980) 


1 2 GM 

E = -v r - 

2 r 



(3) 


where E is the energy, v r is the radial velocity, G is the grav¬ 
itational constant, M is the enclosed mass, r is the radius 
of the shell, and A is the cosmological constant. Adding in 
some tangential motions and making the assumption that 
they are roughly equal to the radial motions, we arrive at 


E = Vr 


GM 

r 



(4) 


Following the prescription of D06, equation (4) is trans¬ 
formed into dimensionless variables, defined as: 


r = 


a y /3 

3GM) r '- 


(5) 


t = 



( 6 ) 


V r 


f^— 

\G 2 M 2 A 


1/6 

Vr 


(7) 


The energy equation is then rearranged to solve for the di¬ 
mensionless time, making use of the fact that v r = dr/dt, 
giving 

t 0 = V 2 f ° ( - L -) V2 dr. (8) 

Jo Vr 3 + 2Er + 2/ 

The difference between equation (8) and the similar equa¬ 
tion from D06 is a multiplying factor of \f2. This inte¬ 
gral can be evaluated analytically for the critical energy 
(E = —3/2, see D06). After performing the integration, it 




Figure 1. Radial and tangential velocities of particles in Object 
1 from Simulation 2 as a function of radius in the top and bot¬ 
tom panels, respectively. The magnitudes of the radial velocities 
are plotted. The vertical dashed lines in both panels mark 7*200 
(density inside is 200p c ,o) for this object. The solid lines show 
average velocities in 100 1 kpc bins. Note the significant tan¬ 

gential motions of particles well outside of the virialized region. 


is useful to define a new variable x(r cs ), as done in D06, 
which will have a different exponent than they found, 

y/3/2 


x(r cs ) = 


1 -f- 2r c s + \J 3r cs (res -1- 2) 
1 + 2r C s — \/3 r cs (r c s + 2) 

x (l + r C s + \Z^cs(^cs ! 2 
Using the fact that 

O — A . i2 ( 3t 

Da. = = tanh 


>) 


-3yf2 


(9) 


( 10 ) 


this variable is then related to Da via 

^ A/1 V I - I ^ 

D a (res) = 


x(r C s) - 1 


x(fcs) + 1 


( 11 ) 


This expression was numerically evaluated by incrementally 
increasing r cs , starting from zero, until DA(r cs ) = 0.7. This 
occured for r cs = 0.75. Framing this in terms of a critical 
parameter, as in equation (1), yields £pi 4 = 1.89 in excellent 
agreement with £ s i m . This is hereafter referred to as the P14 
model. 

Equation (8) was also numerically integrated for ener¬ 
gies other than the critical energy. An energy value was 
selected, then the right hand side was numerically inte¬ 
grated until it equaled the value of the current dimension¬ 
less time, to- This procedure was performed for both the 
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Figure 2. A comparison of the D06, based on the pure SCM, and 
the P14 models. The critical and turn around radii of each model 
are labeled, along with the point at which the two models have 
equivalent results. Noting the locations of the turn around radii, 
at r eq the D06 model would have that shell collapsing, while the 
P14 model would have it still expanding. 


P14 and D06 models, starting with E = 0, incrementing 
until a value of energy where the current radius would be 
zero. The results of these integrations are shown in Fig. 
2. There are several noteworthy features from this figure. 
The turn around radius of the D06 model (rt a ,D 06 = 0.73) 
is only slightly smaller than the critical radius of the P14 
model (r C a,pi 4 = 0.75), and is significantly larger than the 
turn around radius of the P14 model (?t a ,pi 4 = 0.61). There 
is also a point at which the two models cross (r eq = 0.66) 
which lies between the two turn-around radii, meaning the 
P14 model would have that shell expanding, while that of 
D06 would have it collapsing. These features can potentially 
be used to test the performance of the models when com¬ 
pared to simulation data. 

In order to carry out a comparison, the densities asso¬ 
ciated with these shells need to be known, so that they can 
be located in simulation data. The density parameter for 
any shell is given by (D06, equation (17)) 


n 8 



2n A 

ys ' 


( 12 ) 


Since the dimensionless radii of all the shells of interest 
are known, the calculation of the density parameters is 
trivial. The density needed to locate the limits of bound 
structure in the future can also be located easily by us¬ 
ing fl cs (a = 100) = 2 (D06), and multiplying by the future 
(a = 100) value of the critical density. The value for r cs ,B 03 
is determined from the parameter /3 of B03 defined as 


0 


2GM 


(13) 


While they name this the strength parameter, it is actu¬ 
ally just the shell density parameter, which we can see by 


Table 2. Densities associated with various radii from the B03, 
P14, and D06 models. Column 1 labels which radius is being con¬ 
sidered. Column 2 gives the value of that radius in dimensionless 
units. Column 3 gives the shell density parameter. Lastly, col¬ 
umn 4 gives the shell density. The units of density are chosen to 
match the GADGET2 simulation units. 


Radius 

r 

n s 

Ps 

(10 10 h 2 Mq kpc 3 ) 

r cs (a = 100) 

i 

2 

3.88 x 10~ 8 

^cs,D06 

0.84 

2.36 

6.55 x 10~ 8 

^cs,P14 

0.75 

3.32 

9.21 x 10“ 8 

r t a ,D06 

0.73 

3.60 

9.99 x 10" 8 

^eq 

0.66 

4.87 

1.35 x 10“ 7 

r cs,B03 

0.64 

5.28 

1.48 x 10 -7 

rta,P14 

0.61 

6.17 

1.71 x 10“ 7 


rearranging equation (13) as 


3 M 8t tG 
4tt rg 3 H§ ' 


(14) 


The first part of this equation is simply the enclosed mass 
divided by the volume, or the mean mass density, and the 
second part is the inverse of the critical density. From their 
analysis, B03 determine that in order for a shell to have a 
turn around, ft must satisfy 


27 

P 3 ^ — (n m ,o + /3) 2 Ha,o- 


(15) 


Evaluation equation (15) at the equality gives the critical 
density parameter for their model, allowing for the determi¬ 
nation of the density and associated dimensionless radius. 
The results are shown in Table 2. 


4 COMPARISONS WITH GADGET2 
SIMULATION 

In order to rigorously compare all three models, two simula¬ 
tions were performed with GADGET2 (Springel 2005). Simu¬ 
lation 1 was designed to reproduce that performed by D06: 
128 3 particles of mass 3.97 x 10 10 h~ 1 Mg in a 100 h -1 Mpc 
periodic box, fl m ,o = 0.3, S2a,o = 0.7, as = 1, and the soft¬ 
ening length was set to 15 h _1 kpc in co-moving coordinates 
until the present time and then capped at that value in 
physical units. Simulation 2 was larger: 256 3 particles of 
mass 13.4 x 10 1o /i _1 Mq in a 300 /i _1 Mpc periodic box, 
0 = 0.3, fi A ,o = 0.7, as = 1, and the softening length was 
set to 50/i -1 kpc in co-moving coordinates until the present 
time, after which it was capped at that value in physical 
units. Both simulations were evolved from a = 0.02 (z = 49) 
to a = 100. Initial conditions were generated using the N- 
Genic software, also developed by Volker Springel. The pur¬ 
pose of Simulation 1 was to check the methods used here 
for consistency with those of D06. However, given the rela¬ 
tively small box size of that simulation combined with the 
periodic boundary conditions, the long range gravitational 
effects may have been somewhat suppressed. These forces 
might be very influential in terms of the effects of external 
attractors on potentially gravitationally bound structures. 
The larger box size of Simulation 2 should have allowed for 
more realistic long range gravitational forces, better mod¬ 
elling the effects of external attractors. 
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Both simulations were analyzed in the same way. Snap¬ 
shots were taken at the present (a = 1) and in the far fu¬ 
ture when structure formation should decrease significantly 
(a = 100). Structures were identified using the new software 
COLDGaS 1 , Compute Unified Device Architecture (CUDA) 
Object Location Determination for GADGET2 Snapshots, 
developed by the author. As the title suggests, the code uses 
NVidia’s CUDA language to leverage the massively parallel 
capabilities of graphics processor units (GPUs) to signifi¬ 
cantly speed up certain portions. The code first performs 
a grid based density calculation in which particles are as¬ 
signed to cells by dividing each coordinate of the particles 
position by the length of one side of a cell, and rounding 
the result down. This gives three indices corresponding to 
a unique cell. In this manner, the density calculation will 
always have a complexity of order 41V, regardless of the size 
of the grid. The particles are sorted according to which cell 
they belong in, so that only information about occupied 
cells is stored, reducing memory costs for small grid sizes. 
Density peaks are then used to find galaxy cluster like ob¬ 
jects, the most massive of which go on to become the cen¬ 
tres of structures to be studied as potentially gravitationally 
bound. After some further improvements are made, a future 
paper will provide a much more detailed description of the 
software. 

Both the present (a = 1) and future (a = 100) snap¬ 
shots were analyzed with COLDGaS. The most massive fu¬ 
ture objects were matched to their progenitors at present 
based on their positions in the simulation volume. Only ob¬ 
jects that were sufficiently far away from the edges so that 
their radii would not overlap the box boundary were cho¬ 
sen for simplicity. Once objects were matched between the 
present and future snapshots, particle IDs were compared 
to determine which particles remained part of the structure 
into the far future, implying that they are gravitationally 
bound. 


4.1 Simulation 1 

The primary reason for Simulation 1 was to verify that the 
methods used here were consistent with those used by D06. 
In order to check this, the output of that simulation was 
analyzed in the same way as the data of D06. Using the 
above described procedure, 20 objects were selected from 
this simulation. Their future masses ranged between 2.33 
and 19.0 x 10 14 hr 1 Mq. The range in masses at present is 
dependent upon the chosen analytical model, but is in rough 
agreement with the above limits. 

Table 3 shows a summary of the results from analyzing 
the data in the same manner as was done in D06, with the 
results on a per object basis available as online supplemen¬ 
tal material. The four quantities presented are the same as 
calculated by D06 in which the objects are examined on a 
per particle basis. The first two rows of Table 3 show the re¬ 
sults obtained by D06 when they analyzed their simulation. 
Comparing those first two rows with the results obtained 
here in applying the models of B03 and D06 to Simulation 
1, we can see that they agree within the uncertainties. This 
implies that the methods of analysis used here are consistent 

1 For this work, COLDGaS-0.4.3 was used and is publicly available 
at http://davidwpearson.com/some-of-my-code/ 


Table 3. Summary of the analysis of Simulation 1 in the manner 
of D06. Column 1 identifies the model. Column 2 is the ratio 
of particles identified as bound at present that escape to the 
total number selected as bound at present. Column 3 is the ratio 
of particles identified as bound at present that remain part of 
the structure to the total number selected as bound at present. 
Column 4 is the ratio of particles that are part of the structure 
in the future but were not selected as bound at present to the 
total number selected as bound at present. Lastly, Column 5 is 
the ratio of mass that is bound in the future to the mass selected 
as bound in the present. 


Model 

A 

Per cent 

B 

Per cent 

c 

Per cent 

D 

Per cent 

B03“ 

9.9 ±4.2 

90.1 ±4.2 

13.3 ± 12.4 

103.4 ± 14.3 

D06 6 

28.2 ± 13.0 

71.8 ± 13.0 

0.26 ±0.23 

72.0 ± 13.1 

B03 

8.7 ±4.0 

91.3 ±4.0 

20.0 ±22.0 

111.3 ± 19.4 

D06 

19.5 ± 7.0 

80.5 ±7.0 

1.6 ±2.5 

82.0 ±7.0 

P14 

13.7 ±6.0 

86.3 ±6.0 

6.9 ± 10.2 

93.1 ± 10.1 


“Data in this row is from Table 2 of D06. 
■'Dal.a in this row is from Table 1 of D06. 


with those used by D06. As expected, the results of the P14 
model fall between the B03 and D06 models. Both the B03 
and P14 models do well at finding the final bound mass of 
the objects, while the model of D06 predicts far more mass 
being bound. However, all of the model’s results agree with 
each other at the 2 a level, indicating that this is not the 
most informative way of examining their performance. 

For all of the analytical models, the dimensionless crit¬ 
ical radius should tell us the relative size of the present 
critical shell compared to its final critical radius, since all 
models should approach a dimensionless radius of r = 1 as 
t —> oo. For this reason a comparison was made between the 
present and future critical radii, to see if any one model per¬ 
formed better. For the model of B03 the critically bound 
shell should be about 64 per cent of its final size, and in 
Simulation 1 it is on average 62.6 ± 3.4 per cent of its final 
size. For the P14 model, the critically bound shell should 
be about 75 per cent of its final size, and in Simulation 1 it 
is on average 77.3 ± 2.8 per cent of its final size. Lastly, the 
model of D06 predicts the critically bound shell to be about 
84 per cent of its final size, while in Simulation 1 it is found 
to be, on average, 90.2 ± 2.9 per cent of its final size. Tak¬ 
ing only the central values into account would suggest that 
the model of B03 is slightly conservative, the P14 model is 
slightly optimistic, and the D06 model is much more opti¬ 
mistic. However, for all of the models, the predicted ratio is 
within 2cr of the average ratio found in Simulation 1, mean¬ 
ing that we cannot learn much from this analysis. 

It is well known that when two structures merge, some 
of the mass of those structures is ejected (White 1978; 
Colpi, Mayer & Governato 1999). Combining this with the 
fact that all models agree with each other at the 2er level, 
strongly suggests that trying to examine structures from 
simulation on a per particle basis does not make much sense. 
What is of more interest is being able to identify which clus¬ 
ters may be gravitationally bound to each other, and which 
are not. This implies that a more useful way of assessing 
the models’ performance would be to visually examine the 
objects, to see which model does the best at selecting the 
clusters that are bound to the central dominant cluster. By 
using the IDs of particles that are part of the object in 
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the future snapshot, it is possible to pull the coordinates 
of those particles from the present snapshot data. Plotting 
these particles will then give a good idea of what is bound 
to the central dominant cluster. Because of mass being lost 
during mergers, the plots containing only the future parti¬ 
cles will seem less dense than the plots from the analytical 
models. 

Fig. 3 shows the particles selected as potentially bound 
by the three analytical models applied purely spherically 
(only from the central cluster), and the particles that are 
part of the structure in the future frame, excluding particles 
ejected due to mergers. Also included are plots of the par¬ 
ticles selected as potentially bound by the models of D06 
and P14 applied in a ‘non-spherical’ manner described be¬ 
low. Comparing the various panels to the panel containing 
the particles that are part of the final structure (panel d), 
we can get a good idea about the performance of the ana¬ 
lytical models. In the purely spherical mode, the model of 
D06 identifies two substantial clusters as bound to the cen¬ 
tral cluster, and the P14 model only includes one of these 
(the cluster’s centre of mass is inside the P14 critical radius 
by ~ 6 / 2 A 1 kpc), while the model of B03 does not include 
either of them. The one just inside the critical radius of 
the P14 model ends up as bound to the structure, while 
the other cluster does not. The last two panels of the fig¬ 
ure show the ‘non-spherical’ application of the D06 and P14 
models. Here the predictions are made by first applying an 
analytical model to a central dominant cluster, and then to 
any other clusters within the central clusters critical radius. 
Note that panel f (the ‘non-spherical’ P14 model) seems to 
most closely approximate what becomes part of the final 
structure. 

Visual examination of the other objects from Simula¬ 
tion 1 fail to yield any additional insights to the models. 
For figures similar to Fig. 3 for all of the objects from this 
simulation, please see the online supplemental material. 

So far, the results have been rather unsurprising. The 
D06 model seems to be optimistic as to the extent of bound 
structure, while the B03 model seems conservative to the 
point of missing parts of the structure that are bound. The 
P14 model falls between the two, and detailed examination 
of the objects hints that it is slightly on the conservative 
side, though it seems significantly less likely to miss bound 
portions of the structure than the B03 model. 

4.2 Simulation 2 

The larger box size of this simulation should ensure that it 
better approximated the effects of long-range gravitational 
forces, which may have been underestimated by Simulation 
1. The analysis of this simulation was undertaken in the 
same manner as for Simulation 1, again selecting 20 ob¬ 
jects for study with masses in the future ranging between 
1.83 and 7.67 x 10 15 /i _1 Mq. As before, the mass range 
at present will vary depending on the choice of analytical 
model, but it is noteworthy that the most massive object 
in this simulation is ~ 10 16 h _1 Mg at present according to 
the P14 and D06 models, similar to the estimated masses 
of the SSC and CSC. 

Table 4 shows the results of analyzing the data as was 
done in D06. From this table, it is clear that the results are 
virtually identical between the two simulations. Again we 
see that the results of applying the B03 and D06 models to 


Table 4. Summary of the analysis of Simulation 2 in the manner 
of D06. Columns are the same as in Table 3. 


Model 

A 

Per cent 

B 

Per cent 

c 

Per cent 

D 

Per cent 

B03“ 

9.9 ± 4.2 

90.1 ±4.2 

13.3 ± 12.4 

103.4 ± 14.3 

D06 6 

28.2 ± 13.0 

71.8 ± 13.0 

0.26 ±0.23 

72.0 ± 13.1 

B03 

9.7 ±2.8 

90.3 ±2.8 

19.3 ± 15.8 

109.6 ± 14.4 

D06 

19.3 ± 3.8 

80.7 ±3.8 

3.5 ± 5.7 

84.2 ±6.9 

P14 

13.9 ±3.1 

86.1 ±3.1 

7.1 ±8.0 

93.2 ±8.2 


“Data in this row is from Table 2 of D06. 
^Data in this row is from Table 1 of D06. 


the simulation performed by D06, agree with the results of 
applying those models to Simulation 2. The results of the 
P14 model fall between the other two, and agree with the re¬ 
sults of that model from Simulation 1. The ratios of present 
radii to final radii are 62.6 ± 2.6, 77.1 ±2.1, and 89.2 ± 2.4 
per cent, for the B03, P14, and D06 models, respectively. 
These again agree with the theoretical expectations at the 
2 <j level. For the statistics tracked by D06 all of the mod¬ 
els agree with each other at the 2<r level. Results on a per 
object basis are provided in an online supplement. 

Performing another visual examination of the objects 
gives the greatest insights into the accuracy of the models. 
Fig. 4 shows Object 3 from Simulation 2. It can be seen that 
the D06 model selects a cluster that is missed by the P14 
model, and that this cluster is part of the structure in the fu¬ 
ture. The missed cluster’s centre of mass lies ~ 1.5 It -1 Mpc 
outside of the critical radius of the P14 model, making this 
the biggest failure of that model. It is worth noting that 
the central cluster is significantly elongated in the general 
direction of this missed cluster, giving it a very ellipsoidal 
shape. Since the density calculation used the mass enclosed 
by spherical shells, the non-spherical shape of the central 
cluster could explain the failure of the P14 model. It is also 
worth noting that the B03 model significantly underesti¬ 
mates the extent of the structure. For similar figures of all 
objects from this simulation, please see the online supple¬ 
mental material. 

The velocity structure of objects in Simulation 2 was 
also examined to look for the location of the turn around 
radius. To do this, the centre of mass velocity for the cen¬ 
tral cluster was found and subtracted from the velocities 
of all the particles associated with that particular object. 
Then the velocities were broken into their radial and tan¬ 
gential components. Fig. 5 shows the radial velocities of par¬ 
ticles around Object 1 in Simulation 2. The solid line repre¬ 
sents the average radial velocities of particles in 100 h _1 kpc 
bins. From left to right, the vertical lines represent the turn 
around radii based on the P14 model, binned averages, and 
the D06 model. The P14 model underestimates the turn 
around radius, while the D06 model overestimates. The 
same trend is observed in the other objects as well, though 
for some, locating the turn around radius from the binned 
average is difficult due to clusters near that radius. On aver¬ 
age, the P14 model underestimates the location of the turn 
around radius by ~ 16 per cent, while the D06 model over¬ 
estimates its location by ~ 7 per cent. This again suggests 
that the P14 model is slightly conservative, while the D06 
model is slightly optimistic. 
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Figure 3. Three-dimensional plots of the particles associated with Object 1 from Simulation 1. a.) All particles pulled from the present 
(a = 1) snapshot based on the D06 cutoff density, b.) All particles pulled from the present snapshot based on the P14 cutoff density, 
c.) All particles pulled from the present snapshot based on the B03 cutoff density, d.) All particles that are part of the structure in the 
future (a = 100) snapshot in their present (a = 1) locations, e.) All particles pulled from the present snapshot based on the D06 cutoff 
density applied to the central cluster and two additional clusters, f.) All particles pulled from the present snapshot based on the P14 
cutoff density applied to the central cluster and the lower of the two additional clusters. Comparing panels d and f shows that the model 
of P14 applied in the ‘non-spherical’ manner does well at locating what becomes part of the final structure (i.e. what is gravitationally 
bound). 
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Figure 4. Three-dimensional plots of the particles associated with Object 3 from Simulation 2. a.) All particles pulled from the present 
(a = 1) snapshot based on the D06 cutoff density, b.) All particles pulled from the present snapshot based on the P14 cutoff density, 
c.) All particles pulled from the present snapshot based on the B03 cutoff density, d.) All particles that are part of the structure in the 
future (a = 100) snapshot in their present locations. 


5 APPLICATIONS TO THE CSC & THE 
VIRGO CLUSTER 


It was recently shown by PB 2 that five rich Abell clusters 
in the CSC are likely part of a gravitationally bound super¬ 
cluster core. This provides an opportunity to test methods 
of applying these analytical models to a real structure in 
our Universe. The Fundamental Plane (Djorgovski & Davis 
1987) distance estimates published by PB 2 , combined with 
the right ascensions and declinations allow for a reliable 
three-dimensional map of the CSC. In principle, it should 
be a straightforward task to calculate the critical radii of 
the various models by solving equation (1), at the equality, 
for r o yielding 


r o 


( 0.7 2 M obj \ 1/3 
U x 1O 12 M 0 ) 


h _1 Mpc 


(16) 


where M a bj is in units of fi _1 Mq. The only problem is 
that if only the virial mass is known, the critical radii will 
be underestimated, as additional mass outside of the virial 
radius will not be included. Rines & Diaferio (2006) exam¬ 
ined the caustics (Diaferio 1999) around numerous clusters 


in the fourth data release of the Sloan Digital Sky Sur¬ 
vey (SDSS) finding that, on average, the mass within turn 
around is 2.19 ± 0.18 times the mass within the virial radius. 
Anderhalden & Diemand (2011), by examining simulation 
snapshots before the present, located all particles that have 
collapsed at present but expanded beyond the virial radius 
after their first pericentric pass, finding that this would in¬ 
crease the mass of a halo by ~ 25 per cent compared to the 
viral mass. Since their work did not include all the mass 
within turn around, this result is consistent with the work 
of Rines & Diaferio (2006). 

Using the objects in Simulation 1, M 200 , the mass in¬ 
side a radius enclosing 200 times the critical density, was 
used to calculate the critical radii via equation (16) and 
compared with the value determined exactly using the en¬ 
closed density of particles in the simulation. On average the 
B03, P14 and D06 models using M 200 predicted radii that 
were 77.9 ± 5.2, 73.6 ± 6.9, and 73.9 ± 7.6 per cent of the 
values determined by the cutoff densities, respectively. Sim¬ 
ilar results are found using the objects in Simulation 2 with 
75.1 ± 6.7, 71.2 ± 7.1, and 72.0 ± 7.2 per cent, respectively. 
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Figure 5. Radial velocity versus distance from the central clus¬ 
ter for Object 1 in Simulation 2. The large dashed vertical line 
marks the turn around radius predicted by the P14 model. The 
dot-dashed vertical line marks the turn around radius predicted 
by the D06 model. The gray dots are the radial velocities of indi¬ 
vidual particles, and the solid line is the average radial velocities 
in 100 h -1 kpc bins. The small dashed vertical line marks where 
the average radial velocities transition from negative to positive, 
approximating the location of the true turn around radius. 

The virial mass alone gives critical radii for A2065 
of 7.25to.63 , 8.45to.74, and 9.89t£ge /i _1 Mpc for the B03, 
P14, and D06 models, respectively. The confidence intervals 
come directly from the uncertainty in the mass estimate of 
A2065 (2.33±osg x 10 15 h _1 Mq). If we assume that these 
are underestimates as indicated by the simulation data, the 
critical radii become 9.3+J'f, 11.4+ 2 'g, and 13. h _1 Mpc 
for the B03, P14, and D06 models, respectively. Here the 
confidence intervals come from the uncertainty of the de¬ 
gree of underestimation. It is worth noting that these re¬ 
sults agree with the findings of Rines & Diaferio (2006). 
If the virial mass of A2065 is increased by a factor of 
2.19 ± 0.18, the critical radii become 9.41 ± 0.26, 11.0 ± 0.3 
and 12.8 ± 0.3 ft -1 Mpc for the B03, P14, and D06 mod¬ 
els, respectively. Looking at the separations that come from 
the FP distances given in Table 5, we find that the model 
of P14 would include A2056, A2061, and A2089 inside the 
critical radius of A2065. However, A2061 and A2067 are 
only separated by 2.37 h -1 Mpc and the mass of A2061 
is 9.9l2 g x 10 14 h -1 Mq giving it a critical radius of ~ 
9 hT 1 Mpc in the P14 model, taking into account the un¬ 
derestimation from the virial mass. This places A2067 well 
inside the critical radius of A2061 and as seen in section 4, 
everything that is bound to A2061 should also be included in 
the structure. This would imply that A2056, A2061, A2065, 
A2067, and A2089 should be a gravitationally bound struc¬ 
ture, which is in complete agreement with the findings of 
PB 2 . The model of D06 would include all five of the above 
clusters inside the critical radius of A2065 alone. 

The above procedure has the advantage of not needing 
any details about the density profile of the structure. How¬ 
ever, one has to assume that using M 200 from the simulation 
data is approximating the virial mass. Given the agreement 
with Rines & Diaferio (2006), this is likely not a bad as¬ 
sumption, but testing it further is still desirable. In apply¬ 
ing the caustics method to a supercluster (Reisenegger et al. 
2000), the mass is determined as a function of radius. PB 2 
applied the caustics method to the CSC, giving mass as 



Figure 6. Density versus radius as determined from an appli¬ 
cation of the caustics method to the CSC (PB 2 ). The solid line 
are the data from the caustics calculation, and the dashed line 
is a power law fit. Since the ordinate values are normalized by 
the present value of the critical density, they are the shell density 
parameters. 

a function of radius which can be simply transformed into 
density as a function of radius. Fig. 6 shows the density nor¬ 
malized by the current value of the critical density (i.e. the 
density parameter) as a function of radius for the CSC from 
the caustics method. This was fit with a power law yielding 

fi s = (596.6 ± 0.2)r (_1 " 293±5xlo_5) . (17) 

This can then be solved for r using the critical shell density 
parameters given in Table 2. Doing this we obtain 10.7, 13.5, 
and 16.1 h~ x Mpc for the B03, P14, and D06 models, respec¬ 
tively. Taking the cluster separations from Table 5, the P14 
and D06 models again predict that A2056, A2061, A2065, 
A2067, and A2089 are part of a gravitationally bound struc¬ 
ture. 

As a further test, the mass of the Virgo cluster 
within turn around was estimated using the P14 model. 
Karachentsev et al. (2014) examined the infall velocities of 
galaxies around the Virgo cluster. From their analysis they 
were able to determine that the turn around radius for Virgo 
is 7.2 ±0.7Mpc (they assumed Hq = 72 kms -1 Mpc -1 ), 
and using the pure SCM (i.e. the D06 model) they estimated 
the mass of the Virgo cluster to be (8.0 ± 2.3) x 10 14 Mq. 
They make note of the work of Rines & Diaferio (2006), 
but instead use the work of Anderhalden & Diemand 
(2011) to justify their findings as being in agreement 
with the virial mass of Virgo, (7.0 ± 0.4) x 10 14 Mq. Since 
Karachentsev et al. (2014) found the turn around radius, 
it would seem more appropriate to use the findings of 
Rines & Diaferio (2006), which would place the mass of 
Virgo within turn around at (1.5 ± 0.2) x 10 15 Mq based 
on the virial mass. Using the turn around density of the 
P14 model, a mass can be estimated for the Virgo cluster 
via 

Mta = 7.17 x 10 12 h 2 Mq. (18) 

With the turn around radius found by Karachentsev et al. 
(2014), the mass estimate for Virgo would be 
(1.4 ±0.4) x 10 15 Mq, in agreement with the work of 
Rines & Diaferio (2006). 
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Table 5. Relative separations of the clusters in the CSC from FP distances. Each column and row corresponds to a different cluster in 
the CSC. Matching one cluster’s column to another’s row gives their relative separation in units of h ~ 1 Mpc. Individual cluster distances 
have uncertainties of ~ 5 per cent (PB 2 ). 



A2056 

A2061 

A2065 

A2067 

A2079 

A2089 

A2092 

CL1529+29 

A2056 

- 

9.08 

4.63 

10.48 

22.60 

11.18 

16.57 

30.83 

A2061 

9.08 

- 

11.11 

2.37 

24.85 

13.51 

12.89 

28.96 

A2065 

4.63 

11.11 

- 

12.41 

25.09 

8.83 

17.90 

27.55 

A2067 

10.48 

2.37 

12.41 

- 

23.80 

13.27 

10.77 

29.90 

A2079 

22.60 

24.85 

25.09 

23.80 

- 

22.60 

17.51 

50.51 

A2089 

11.18 

13.51 

8.83 

13.27 

22.60 

- 

13.22 

28.63 

A2092 

16.57 

12.89 

17.90 

10.77 

17.51 

13.22 

- 

36.36 

CL1529+29 

30.83 

28.96 

27.55 

29.90 

50.51 

28.63 

36.36 

- 


6 DISCUSSION 

Identifying the extent of gravitationally bound structure 
could be very useful in placing constraints on models of our 
universe (Sheth & Diaferio 2011; Park et al. 2012). Given 
the lack of spherical symmetry, effects of external attrac¬ 
tors and tangential motions due to past gravitational inter¬ 
actions, setting the limits of bound structure is clearly not 
a straightforward task. Here, building on the work of D06, 
a simple modification of the SCM was made, yielding a new 
analytical model (P14) for the extent of bound structure 
in a universe containing cosmological constant dark energy. 
The P14 model was found to agree remarkably well with a 
model informed by the simulations of P13. 

The two previous analytical models for the extent of 
bound structure, B03 and D06, along with the P14 model 
were tested against two GADGET2 simulations, which al¬ 
lowed for objects of a variety of sizes to be examined. By 
visually comparing the structures predicted as bound by the 
various models at present (a = 1) to the structure as defined 
by the particles which remain in the far future (a = 100), 
the B03 model was found to consistently miss portions of the 
structure. The D06 model fairly consistently includes clus¬ 
ters which are not part of the structure in the far future. 
The P14 model performed well in many cases, though it 
does still seem to slightly underestimate the extent of bound 
structure. The fact that the P14 model performed well, and 
that it agrees so well with the P13 informed model suggest 
that those simulations, while not very sophisticated, were 
still able to reliably locate bound structure as long as the 
peculiar motions along the line-of-sight are not significant, 
as they were in the case of the CSC (PB 2 ). This would seem 
to suggest that by giving the clusters a Hubble expansion 
velocity at present, combined with open vacuum boundary 
conditions, the simulations of P13 were in some ways mod¬ 
elling the effects of external attractors. 

In order to truly locate all structure that is bound, it 
is necessary to go beyond the simple spherical models. As 
seen here, after locating all the clusters that are predicted 
to be bound to the central cluster, if one then locates what 
is predicted to be bound to those clusters, it is possible to 
much more accurately determine what is part of the struc¬ 
ture in the far future. This ‘non-spherical’ approach is still 
fairly simple, and straightforward to apply. Future versions 
of COLDGaS will do this automatically when locating po¬ 
tential gravitationally bound objects. Also, Object 3 from 
Simulation 2 (see Fig. 4) seems to suggest that when the 
central cluster itself is not very spherical, it may be neces¬ 
sary to alter the shape of the shells used in the calculation 


of the density to more closely resemble the central cluster. 
Another possible approach may be to build an analytical 
model for the limits of bound structure from an ellipsoidal 
collapse model (ECM). Recently, Rossi, Sheth & Tormen 
(2011) used an ECM along with excursion set theory to 
study the shapes of dark matter haloes. The analytical 
model they derive was shown to more accurately model 
non-linear structure evolution than the Zeldovich approx¬ 
imation. Given the lack of agreement about how to define 
non-spherical gravitationally bound structures from simu¬ 
lation, they instead focus on comparing the axial ratios of 
structures in simulation to those predicted by their model. 
Strictly speaking, their work does not attempt to define the 
limits of bound structure, yet extending such an ECM could 
potentially provide a much more accurate analytical model 
which could inform the COLDGaS software. It is hoped that 
further study of bound structures from simulation, as well 
as incorporating the lessons learned in this work, will lead 
to a version of COLDGaS capable of reliably locating gravi¬ 
tationally bound structures in simulation. 

Three different methods of applying these analytical 
models to real structures in our Universe have been shown. 
The first uses only the cluster virial masses, and makes some 
adjustment based on data from the simulations. In applying 
this method to the CSC, the P14 and DOG models both agree 
with the results of PB 2 . This method is definitely not ideal 
because of the assumption that M 200 = M v i r i a i. However, 
given the agreement of the critical radii adjusted based on 
that assumption, and the critical radii from increasing the 
mass according to the results of Rines & Diaferio (2006), 
this should be a viable method for locating the extent of 
bound structure when very little data exists for the inter¬ 
cluster regions. 

For structures that have been well studied, and have 
ample data available, knowledge gained from a caustics 
analysis may provide a better estimate of the extent of 
bound structure. The caustics analysis of the CSC (PB 2 ) 
was used to yield an estimate of the density profile of the 
structure. With this profile, it was a simple matter of lo¬ 
cating the radii at which the various cutoff densities are 
reached. Using this method, the P14 and D06 models again 
agree with the results of PB 2 . 

Since the turn around densities of the P14 and D06 
models are known, if the turn around radius of a real struc¬ 
ture can be reliably located, it is straightforward to then 
estimate the mass within turn around. Karachentsev et al. 
(2014) did just that with the pure SCM (D06 model). Us¬ 
ing the P14 model, the mass within the turn around radius 
of the Virgo cluster was found to be (1.4 ± 0.4) x 10 15 Mg 
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which agrees very well with the findings of Rines & Diaferio 
(2006), that the total mass of a cluster inside the turn 
around radius is 2.19 ± 0.18 times the virial mass. 

Future work will examine these objects in greater de¬ 
tail, testing deviations from spherical shells for density cal¬ 
culations, performing an in depth examination of their ve¬ 
locity structure, as well as performing mock observations 
to gain a better sense of how these bound structures would 
appear in redshift survey data. Additionally, given that the 
true extent of bound structure seems to lie somewhere be¬ 
tween that predicted by the P14 model and that predicted 
by the D06 model, it should be possible to develop a model 
informed by the GADGET2 simulations performed here. In 
the end, this future work will seek to develop a method 
in the vein of Diinner et al. (2007), similar to the caustics 
method, in order to provide a mass estimate and dynamical 
description of structures. 
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